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1 Introduction 

A number of hematological diseases are characterised by oscillations in the circulating 
density of various types of blood cells. These include chronic myelogenous leukaemia 
(CML), cyclical neutropenia (CN), polycythemia vera (PV) and aplastic anaemia 
(AA). Examples of blood cell counts for CML and CN are shown in figures ^ and |21 

A review of the clinical data, and a discussion of possible mechanisms for the 
oscillations, is given by Haurie et al. (1998). These mechanisms focus on the role of 
negative feedback control on proliferation and differentiation of blood cells within the 
bone marrow, together with time delays due to cell cycling and maturation. There 
are consequently a number of different ways in which oscillations can occur, and one 
object of mathematical modelling of blood cell development is to understand which 
of these effects may be responsible for the oscillations which are seen. 

Blood cells are produced through a process of differentiation from primitive stem 
cells in the bone marrow. These pluripotential stem cells begin to develop along one 
of several different cell lineages, forming blast cells which eventually develop through 
a number of different stages to form the various different kinds of blood cells. The 
most numerous are the red blood cells, or erythrocytes, whose normal density in the 
blood is about 5 x 10 6 cells /d _1 . Their primary function is in transporting oxygen 
to the tissues. Platelets are formed by the fragmentation of megakaryocytes, which 
develop in the bone marrow. Platelets are present at levels of 5 x 10 5 cells /il _1 , and 
their function is in blood clotting. Finally, there are a number of different white blood 
cells, the most common of which are neutrophils (5000 cells /il" 1 ) and lymphocytes 
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Figure 1: Oscillations in white blood cell, platelet and reticulocyte numbers in a 
patient with chronic myelogenous leukaemia. The units are white blood cells, 10 5 cells 
/il _1 ; platelets, 10 5 cells /d _1 ; reticulocytes, 10 4 cells /d -1 . Redrawn from Chikkappa 
et al. (1976). 



(2000 cells /d -1 ), which form constituent parts of the immune system. The normal 
levels of these cells are controlled by a number of mechanisms, and excess or deficit 
of the various cell types defines certain kinds of disease; for example, anaemia refers 
to a low red blood cell count, below 4 x 10 6 cells /d -1 . 

There are a number of features in figures Q and 121 which are of note. In CML, 
there are regular oscillations in white blood cell counts with a long period ranging 
from 40 to 80 days (Fortin and Mackey (1999)). The other cell lines (platelets and 
reticulocytes, i. e. erythrocyte precursors) also oscillate in a similar fashion (figure ^ 
does not show this well: see Fortin and Mackey (1999) for other examples). 

A similar observation is true of cyclical neutropenia. Oscillation periods are of 
order 20 days, during which there is a marked collapse of the neutrophil count to 
vanishingly low levels (Dale and Hammond 1988, Guerry et al. 1973). Other cell types 
oscillate, but only the neutrophils appear to oscillate fairly regularly: oscillations in 
other cell types (e.g., red blood cells, platelets, reticulocytes and lymphocytes) are 
marked by irregularity and high frequency 'noise' (Guerry et al. 1973). This latter 
feature is well illustrated in figure 121 

The purpose of the present paper is throw some light on these observations by the 
study of a model of blood cell proliferation and differentiation. This model is sim- 
ilar to those of previous authors, particularly that of Mackey and Rudnicki (1994), 
and describes the stem cell and developing (blast) cell populations as functions of 
time, age (time through the proliferative cell cycle), and maturation (stage in the 
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Figure 2: Oscillations in neutrophils, platelet and reticulocyte numbers in a patient 
with cyclical neutropenia. The units are neutrophils, 10 3 cells x ; platelets, 10 5 
cells /il -1 ; reticulocytes, 10 4 cells /A' 1 . Redrawn from Haurie et al. (1998). 



differentiation process). Fokas et al. (1991) describe a model with discrete genera- 
tions in the development of blast cells, while Mackey and Rudnicki (1994) develop a 
corresponding continuous model (i.e., developmental stage is a continuous variable). 

In this paper we use a continuous model to describe the development of a single cell 
lineage following the committal of stem cells. Three separate controls are implemented 
in the model, namely the proliferative control of stem cells, the proliferative control of 
developing blast cells, and the peripheral control of stem cell committal by circulating 
blood cell density. We show that variation of parameters in all three control systems 
can cause oscillations, and that the characters of these oscillations are very different. 
This allows us some potential insight into the mechanisms that may be operative in 
some of these dynamic blood diseases. 

2 A model of maturation of blood cell production 

The basic model is similar to that described by Mackey and Rudnicki (1994). It has 
been analysed in various versions by Rey and Mackey (1993), Crabb et al. (1996) and 
Mackey and Rudnicki (1999). A particular feature of these models was the assumption 
of zero maturation rate at maturation state zero. In our formulation of the model, 
we do not make this assumption. 

We consider all cell lineages to consist of populations of two types, proliferative 
and resting phase. These are denoted p and n, respectively, and are functions of age 
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a (time since the inception of the proliferative cell cycle) and maturation m (degree 
of maturation, measured in maturation units (mat), which could be, for example, 
cell generation number). Also, p and n are functions of time t. Thus, we have 
p = p(t, m, a) and n = n(t, m, a). The dimensions of p and n are cells age -1 mat -1 . 

In the cell population, there will be a finite number which are primitive and have 
not begun differentiation. These can not be characterised in terms of p and n at 
m — 0, since the latter are cell densities with respect to a and m. The primitive stem 
cells are characterised by cell densities po(t, a) and n (t, a), such that p da and n da 
are the numbers of primitive stem cells (with m = 0) of age in (a, a + da). 




m (maturation) 



Figure 3: Schematic evolution of cells. Each cell ages as it goes through its cell cycle, 
before dividing and entering a resting (Go) phase; at the same time, the cells mature. 
The time-like variables a (age) and m (maturation) are independent. 



The evolution of the system is illustrated schematically in figure El We suppose 
that cell mortality occurs at a rate 7 (for proliferating cells only), and that cell 
maturation occurs continuously at a rate V (for both proliferative and resting phases). 
We suppose that both 7 and V may depend on maturation stage m, but not on t. 
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Conservation of proliferative cells then implies 



dp + dp + d(Vp) = ^ ^ 

dt da dm 



where the units of V are mat d~ x (maturation units per day). We suppose ()2.1|) 
applies during a cycle of length r (which might depend on m), thus for < a < r; 
then for a > r, the cells in the resting phase satisfy the equation 

dn dn d{Vn) 
dt da dm 

which differs from (12. lj) by the rate of recruitment R back to the proliferative phase; 
resting cell mortality is taken to be zero. Equation ()2.2|) applies for a > r. 

At the end of the cell cycle, a = r, we apply a boundary condition describing the 
conversion of p to n. Mackey and Rudnicki (1994) allow a very general condition, on 
the basis that cells at maturation M can divide to form cells at maturation g(M) < M. 
Specifically, this implies 2p[t, M, r(M)] dM = n[t, g(M),r{g(M)}} dg. If we write 
m = g(M), M = h(m) (so h = g~ r ), then this becomes 

n[t,m,T(m)] = 2p[t , h(m) , r {h(m)}]h' (m) , (2-3) 

where g(m) < m implies h(m) > m. 

A boundary condition for p at a = follows from the recruitment condition (the 
renewal equation) 

p(t,m,Q) = RN(t,m), (2.4) 
where we introduce the total resting cell density 



CO 



N = J nda. (2.5) 

Now we integrate (J2.2|) from a = r to a = oo, taking n — * as a — > oo (which 
is necessary if there are a finite number of cells). We suppose that V = V(m) is 
independent of a and t, and and R = R(t, m) is independent of a. Then 

dN d(VN) 

~dt + = ~ RN + 2p[t ' ^ m) ' T {Km)}}h\m) 1 (2.6) 

adopting ([2.3)1 . 

We need to solve (|2.1|) for p. We use the method of characteristics, and begin 
by applying the recruitment condition (|2.4jl . Specifically, we apply the parametric 
conditions 

t = s, m — [A, a = 0, p = R(s, p,)N(s, //), (2.7) 
valid for s, fi > 0; then the characteristic solution is 

t-S, I y^=t~S, 



p = R(s, fM)N(s, fi) exp 



J s 



(2- 



5 



Define a function v(m,a) by 



dp 



(2.9) 



Then a = t — s, fi = u(m, a). Also dt = dm/V(m) on a characteristic, thus for t > a 
(and also v > 0), 



p(£, m, a) = R[t—a, 1/(171, a)]N[t— a, u(m, a)] exp 
simplifying and putting a = r, we have 

p(t, m, t) = R[t—T, z/(m, r)] N[t— r, z/(m, r)] exp 
for t > r and z/ > 0. Finally, (|2.6|) becomes 



(2-10) 



u{m,T) V(p) 



V[v(m, r)] 
V(m) 



(2.11) 



<9iV 9 

+2//(m) J R[t - r, z/{/i(m), r}]iV[t - r, u{h(m), r}) x 



exp 



Note that 



u{h(m),r} V(p) 

dp 



V[u{h(m),r}} 
V[h(m)\ 



T. 



lu(m,T) V(p) 

It is convenient to define a modified maturation variable £ by 

dp 



£= / 

Jo 



v(pY 



(2.12) 
(2.13) 

(2.14) 



£ has units of time, and indeed it is equal to the elapsed time during maturation. 
Note that v > if £ > r. The lower limit can be chosen for convenience, and allows 
us to fix £ at some reference point; here we choose this to be the initial maturation 
stage (note that this cannot be done if ^(0) = 0). Define also 



h (m) dp 



V(p) 



(note 77 > £ since h > m). Now if 



F(m) = /(£), 



(2.15) 



(2.16) 



then we find 



F[h(m)\ 
F[u{h(m),r}} 



f(v-r). 



(2.17) 
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We change variable from m to £, and define 



v(0 = V{m\ 
M = NV 



(2.18) 



(note that Md^ = Ndm, so that M is cell density in terms of the variable £; the units 
of M are cells d _1 ). After a little manipulation, we find 



dM dM 



-RM + Q, 



where 



Q = 2r)'(£)R[t -T,ij-T]M[t-T,rj- r] exp - f 7 d£ 



(2.19) 



(2.20) 



where we write 7, R and M as dependent on £ rather than m. This equation applies 
if t > r and 77 > r. 

In order to find the form of the source term for t < r or 77 < r, we must solve the 
equation ()2.1j) for p using the initial data from m = and t = 0. If, specifically, we 
have an initial condition 

p = pi(m,a) at * = 0, (2.21) 
then after some algebra we find that 



Q = 2r/ (OpjW - *, t - t v(77 - t) exp - / 7 d£ 



t < r, 77 > i. 



(2.22) 



The definition of Q in t > 77 and 77 < r requires consideration of the stem cell 
evolution, and we now turn to this. Conservation laws for the stem cell densities po 
and n are 



dpo | dp 

dt da 
dn dn 



Of 



da 



-(70 + V )po, 
-(V + R )n , 



(2.23) 



where Vo is the rate of loss of stem cells to maturation, Rq is the stem cell recruit- 
ment rate from the resting phase, and 70 is the mortality rate of stem cells in the 
proliferative phase. We allow _R , V and 70 to depend on t, but we assume they are 
independent of a. Note that V 7^ 0, indeed V 7^ V(0), as the units of V and V are 
not even the same: V has units of mat d , while Vo has units of d _1 . Note also that 
Po and n have units of cells age -1 (unlike p and n). 

The primitive loss to maturation must balance the source for p and n at m = 0, 
thus 

V0P0 = (Vp)\ m=0 , V n = (Vn)\ m=0 , (2.24) 

and the units are consistent. 

Analogously to (|2.4|) and (j2.Hjl . we have 



Po(t,0) 
n {t,r) 



Ro(t)N {t), 
2po{t,r), 



(2.25) 
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where 

Integration over a now yields 
and 



N = n da. 



-VqNq - R N + 2 Po 



(NV)\ m=0 = N V . 
In order to find po we must solve 

^ + ^ = -(70 + ^0, 
with parametric initial conditions: 

Po=p 00 (a), a = a > 0, t = 0, 
Po = R (s)N (s), a = 0, t = s > 0. 



(2.26) 

(2.27) 
(2.28) 

(2.29) 



For t > a, the solution is 

Po = Ro(t — a)N (t — a) exp 
whereas for £ < a, 



t—a 



[7o(0 + Vb(*')]d*' 



= Poo (a - t)exp 



(2.30) 



(2.31) 



(2.32) 



Putting a = r, we find 
dN 



dt 



(R + Vo)N + 2R (t - T )N (t - r) exp 



t-T 



ho(t') + V (t')}dt> 



t > T, 

(2.33) 

which prescribes the control system for iVo, analogously to that of Mackey (1978). For 
t < t, the equation for N involves the initial condition for p , and we can equivalently 
simply prescribe initial data for iVo there. 

Finally, we complete the definition of Q in (|2.19J) by solving (|2.1|) using the initial 
data on m = 0: 



m = 0, a = a>0, t = s > 0, V(0)p = V {s)p (s, a). (2.34) 



We find 



Pit^a)^-^-^- ^ 



from which it follows that 



Q = 2rf(0V Q (t - v)Po[t ~V,r-v] exp [- /" 7 d£ 

[ Jo 



e<a, £<t, (2.35) 



t > T), T] <T. (2.36) 
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Along with ()2.2(Jj) and ()2.22|) . this completes the definition of Q for t > 0, 77 > (and 
thus £ > 0). In summary, 



Q 



2ri'(£)R[t -t,t]- r]M[t -t,T]-t] exp 
2r]'(^)pj[r] — t,r — t]v{r] — t) exp 



Tj — T 



T)-t 



7^ 



t > T, T] > T, 



t < T, T] > t, 



I" rv 

V(0 W - t/W* - 77, ^ - 77] exp -/ jd£ 
*. L Jo 



t > 77, 



7] < T. 

(2.37) 

The two equations (|2.33|) and ()2.19|) are coupled through ([2.28)1 . which provides 
the requisite boundary condition for M at £ = 0: 



M = V N at £ = 0. 



(2.38) 



We see that by an appropriate consideration of the primitive stem cells, we derive a 
coherent model which does not require V(0) = 0. 

Many authors (for example, see Rey and Mackey (1993) and Dyson et al. (1996)) 
study the differential equation ()2.12|) for N, under the assumption that V does tend 
to zero as m — > 0, for example, 

V = rm. (2.39) 

The reasoning behind this is that, if primitive stem cells mature at a finite rate, then 
such cells will be immediately lost to m > 0, which makes no physiological sense, 
because the cell population then inexorably disappears. Only by choosing V(0) = 
can we allow primitive stem cells to endure. In the present version of the model, it is 
also possible to allow V(0) = 0, for example, ()2.39|) would then imply 



M -> as i 



-00. 



(2.40) 



The sensitivity of the solution to this condition has led to the idea that the system 
may have unstable and even chaotic solutions (e.g., Crabb et al. (1996)), because of 
the degeneracy of the equation at m — 0. Our considerations here suggest that the 
requirement that V(0) = is inaccurate, because it does not properly address the 
biological question of how the primitive stem cells should be described. 



3 Dimensionless model 

How we solve the model depends on the complexity of our assumptions about 7, R, 
r) and Rq. In the remainder of this paper we will assume g(m) = m, thus 77 = £, and 
that 7 and 70 are constant. The equation for the maturing cells, (|2.19|l . is 



dM dM „ . . 

_ + _ = _ BM + Q , p.!) 

and is a hyperbolic delay-partial differential equation. Figure H] shows the regions 
where the different definitions of Q apply. In regions II and segment (a) of region 
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a 



b 



III 



II 



T| = T 



t = T 



Figure 4: Regions of different definitions of Q in (|2.37p . Regions I, II and III corre- 
spond to the first, second and third definitions of Q and their locations of validity in 
the (t, 77) plane. The vertical dotted line in region III (where Q is defined in terms 
of po) separates the region (a) where (I2.32j) applies (to the left) from that (b) where 
(I2.31J1 applies (to the right). 



Ill, that is for t < r and all rj = £ > 0, Q depends on the initial data, either pj (in 
II) or p 00 (in 111(a)). Thus we may equivalently simply choose instead to prescribe 
M in < t < t, and this we do. In fact, since £ is finite, the part of the solution 
which depends on this initial data will wash out of the system in a finite time. It is 
therefore apparently of little concern. 

We therefore confine ourselves to consideration of the definition of Q in regions I 
and III(o); with the assumptions we have made, we find that for t > r, 



2e-i T R{t - r, £ - r]M[t - r, £ 



Q 



2 e -7o(T-O e -7£ 



exp 



t-T 



V (tf) dtf 



V (t-ORo(t-r)N (t~r) 



(3.2) 

We thus have to solve ()3.1|) with ()3.2|) in t > r, with the boundary condition (|2.38j) 
on £ = 0, and prescription of an initial function for M in < t < r. 

A principal issue of focus is how the recruitment rates R and Rq and the committal 
rate Vq depend on M, Nq and £. There is very little to constrain our choice. In what 
follows, we assume Rq = Rq(Nq) (stem cell proliferation is controlled by stem cell 
density). We follow Mackey and Rudnicki (1994) in supposing that R depends on the 
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total differentiating cell population M, where 

rip 

M(t) 

£f being the time of final maturation, 



mp 



dp 



v( P y 



(3.3) 



(3.4) 



and m = mp at full maturity. We suppose that the rate of committal Vq should 
depend on the peripheral blood cell count, B, thus V = V (B). A simple model for 
B is 

d R 

— = M| fF - lB B, (3.5) 

where 7# is the specific decay rate of the peripheral blood cells, and the source term 
M|g F is the delivery rate to the blood from the maturation phase cells. Peripheral 
control models of similar type have been studied by Bernard et al. (2003). Assump- 
tions of this type are liable to be important in the evolution of diseases such as 
cyclical neutropenia, which is thought to be due to an instability in the peripheral 
control of stem cell committal. In addition, it is likely that other controls affect rate 
of apoptosis, maturation rate, cell cycle time, and so on. 
The equation for Ao (I2.33J) now takes the form 



N = -[V + R (N )}N + 2e-^ T J Ro(Ao T )A OT exp 



t-T 



V [B(t')]dt' 



(3.6) 



where A 0r = N (t — r). This is precisely the model of Mackey (1978) if V is constant, 
and has been studied by Fowler and Mackey (2002) in the limit 



V t < 1, 



(3.7) 



when it is shown that relaxation oscillations will occur for a further parameter p = 
[2e~ < - 70+v ' ' |T — l]/VoT within a certain 0(1) range. (Note that in the notation of 
Fowler and Mackey's model, 7 = 70 + V0, 6 — Vo.) When such oscillations occur, they 
will propagate through the maturing cells; however we show in this paper that the 
resultant amplitude of oscillations of mature blood cells is small unless amplification 
also occurs during maturation. 

We can write (J3.2j) in abbreviated form as 



2e~^ T R{M T )M T , 



Q 



2 e -7o (•?--?) g-7£ 



cxp 



t-T 



V [B(t')]dt' 



V [B(t - 0]Ro(N Qt )N 0t , 



(3i 



where M T = M(t - r), and M T>r = M(t - r, f - r). 

We non-dimensionalise the model by following the analysis of Fowler and Mackey 
(2002), which motivates a choice of scales for the variables as follows: 

/•( - 7. M ~ M*, A,, .Y,;.S. Q~— , 

T 
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Ro = R* h , R = R*h, V = V*v , B 



M* 

IB 



(3.9) 



where Rq, R*, Vq, Nq, M* are determined by the control functions (so that they 
are 0(1) functions of 0(1) variables). For example, Mackey (1978) chooses for Rq 
the Hill function 



Rq 



l + (N /9) n - 

In this case, we would choose N£ = 9 and ho is the Hill function 

ho(S) 1 



1 + S r 



The dimensionless stem cell equation is 



S = b 
where 



(1 + X )h (S l )S 1 exp £ 1- MBit')} dt'\\- h S 



t-i 



e = V*r, A = 2e~^ +v ^ - 1, b = R*r. 
The dimensionless form of ()3.8j) is 

dM DM 



where 



Q 



dt + d£ 



b(i + \)h(M 1 )M 1<1 , e>i, 



—bh(M)M + Q, 



(3.10) 
(3.11) 

-e v S, (3.12) 
(3.13) 

(3.14) 



z/6 (l + A )e- a? exp 



e h- £_*v [B{t)]d1f 



v [B(t - ()]h {3 1 )3 1 , £<1, 
(3.15) 



in which 



N*V* 
M* 



b = R*r, X = 2e~ 7T - 1, a = (7 - j )r. 



(3.16) 



This is analogous to the scaling used by Fowler and Mackey (2002). The boundary 
condition for M is 

M = vv S at £ = 0, 



and if 

then 

where 



M = M f at £ = £/, 



SB = M f - B, 



JbT t 



(3.17) 
(3.18) 
(3.19) 

(3.20) 



This completes the statement of the dimensionless form of the model. 
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Parameter values 



The equation ()3.12|) is exactly that studied by Fowler and Mackey (2002) (if Vo = 1). 
However, their model can also be interpreted as a 'lumped', or compartmentalised, 
version of 1)3.14)1 for the maturing cells. One way of enabling this is if we make the 
special assumption that the maturation rate V — > as both m — > and m ^ m^, as 
also assumed by Mackey and Rudnicki (1994). In this case the range of £ is (— oo, oo), 
and we have M — > at both limits. Then integration of (|3.14|) over £ again leads 
to an equation of the form of (j3.12j) . In the present paper we assume V is finite at 
m = rrip, i.e., the mature blood cells are delivered to the bloodstream at a finite 
rate, and this is then the essential difference between the models with and without 
maturation. 

In estimating the parameters, we follow Fowler and Mackey (2002) in choosing 
t ~ 2 d, and we suppose that proliferative control is effected at typical rates R* ~ 
Rq ~ 2 d -1 . We suppose apoptosis rates are of order 7 ~ 70 ~ 0.2 d -1 , and that 
committal rates are of order V^* ~ 0.05 d _1 , and from these we find 

b~b ~4, A~A ~0.3, e ~0.1. (3.21) 

The parameter a is not independent of the others, as 

a = e + In (j^y) , (3-22) 

and plausibly a ~ e . 

The remaining parameters are 6, v and £/. For 5, we assume a half life (7s 1 ) of 7 
hours, appropriate for neutrophils (but, for example, certainly not for erythrocytes); 
then 

5-0.15. (3.23) 

We can get some sense of the size of the remaining parameters v and £/ by con- 
sidering the nature of stem cells. These are difficult to isolate; indeed it is not yet 
clear whether genuine stem cells have ever really been isolated. The reason for this 
is that there are few of them, and maturing cells will typically undergo about (or at 
least) twenty divisions before emerging as mature blood cells. A typical numerical 
estimate for the total number of blast cells is 10 12 per kg body weight, while for stem 
cells, a corresponding estimate is 10 6 (Bernard et al. (2003), Mackey (2001)). If this 
is the case, then it successively implies that the parameter v in (I3.16|) is very small 
(w 10~ 6 ), and therefore also that the maturation time is long. Typical estimates 
of £f ~ 10-20 days are consistent with values of £/ ~ 5-10, and in fact the small 
parameter ^J 1 then plays the role corresponding to that of the small parameter e in 
Fowler and Mackey's (2002) analysis. 

Steady state 

To elaborate this discussion, we now describe the steady state. For simplicity, we 
ignore the distinct definition of Q in £ < 1, and extend the definition in £ > 1 back 
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to £ = 0. The steady solution of ()3.14|) and ()3.15|) is, with v = S = 1, 



M = ve s Z, (3.24) 
where s is the unique positive solution of the pair 

s = bh(M) [(1 + A)e~ s - l] , 

M = - (e s ^f - l) . (3.25) 

(Uniqueness follows from the fact that M is monotonically increasing with s, hence 
h(M) is monotonically decreasing with s, hence bh(M) [(1 + A)e~ s — 1] is monotoni- 
cally decreasing with s, while evidently s is increasing.) We can see that s < ln(l + A), 
and s will be close to this value if b is large. Note also that by choosing A and e to 
have certain specific values which depend on A, b, v and £/, this solution consistently 
extends back to £ = 0, even allowing for the distinct definition of Q in £ < 1. 

Numerical solutions do confirm the exponential variation of M with £. In general, 
it is found that M decreases for < £ < 1, before subsequently increasing. 



4 Periodic solutions 

We are interested in finding whether periodic solutions can occur. There are three 
different controllers in the model, and thus three different ways in which oscillations 
can occur: these are described below. We define a default reference set of parameters, 
and these are given in table [TJ They are those suggested by independent estimate, 
except that we take v = 10~ 2 rather than 10 -6 . This is partly for numerical expedi- 
ency, as smaller values of v require larger and thus longer computation times, and 
also because the value of v is not well constrained. 



Symbol 


Typical value 


n 


3 


£o 


0.1 


Ao 


0.3 


b 


4 


V 


10~ 2 


b 


4 


A 


0.3 




5 


5 


0.2 


V* 


1 


v' 


1 



Table 1: Default parameter values. 



14 



We choose the Hill function controller (|3.11|) for both h and ho, thus 

h{M) = TVW' h o(S) = Y^, (4.1) 

and we take the peripheral controller vq to have the form 

v = [v* - v'B} + , (4.2) 

with default values of the amplitude and slope parameters to be v* = v' = 1. The 
choice of a threshold in ()4.2j) is motivated by the observation that neutrophil popula- 
tions can dwindle to zero in cyclical neutropenia, which would appear to require zero 
production for sufficiently high blood cell counts. 

With this choice of the controller functions and using the default parameters, the 
steady state is stable. Instabilities arising from parameter variations are described 
below. 



Numerical method 



We have to solve the two ordinary differential equations (|3.12j) and (|3.19j) . and the 
partial differential equation (|3.14j) . The solution is complicated by the presence of 
the integral in (|3.12j) . We define 

rt 



U = Sexp 



eo J v [B(t')]dt' 



(4.3) 



and then S and U satisfy the pair of equations 

U 



s = s 



u 



- £qv 



U = b Q [(1 + X)e a h (S 1 )U 1 - h (S)U] . (4.4) 

On the assumption that S remains bounded, U grows exponentially as U ~ exp(v t), 
where Do is the mean of vq. This is likely to cause difficulty in numerical solutions, and 
this can be reduced by using the algebraically growing function W = In U, whence 



S = S 



W - EqVq 



w 



b Ul + A)/io(5 1 ^ a+w/l " H/ 



(4.5) 



>ije~' •' - h (S) . 

In our numerical solutions, we solve ()4.5j) using the second order accurate improved 
Euler method, and we similarly solve (|3.14jl along the characteristics £ — t = rj, on 
which the function Q takes the form 

b{l + X)h{M 1 )M 1} £>1, 

(4.6) 

6o(l + X)e a ^M ( V )h (S 1 ) exp [W x - W{rj)\ , £ < 1, 

where Mo(r)) = M(r),0) and Mi = M^i, i.e., M delayed by one along the character- 
istic. 

Accurate solutions are obtained with a time step At = A£ = 0.05, and these were 
checked aginst values At = A£ = 0.02 (which are used to give the figures). 



Q 
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1 1 1 1 1 1 

400 410 420 430 440 450 



Figure 5: Default parameter values, except that A = 0.05. Stem cell oscillations are 
induced, without any significant effect on blood cells. 



Stem cell oscillations 

Oscillations in the primitive stem cell population will occur for a finite range of 
the parameter Xq/eq, as described by Fowler and Mackey (2002), when v = 1. 
For the default values of b Q = 4, n = 3, the approximate range of instability is 
0.5£o ^ Xq < 1.5eq, and this is modified in an obvious way when the peripheral con- 
troller alters the value of vq. Figure El shows the oscillations which occur in the stem 
cell population when Ao is reduced to 0.05. It is an interesting fact that these oscil- 
lations are hardly manifested in the blood cell population. The apparent reason for 
this is that the small value of v means that oscillations in M are small, and therefore 
also in Mf, because small perturbations propagate stably down the maturation axis. 
The blood cell population is therefore stable, and B « Mf. 



Proliferation-controlled oscillations 

We use the term proliferation-controlled oscillations to refer to oscillations induced 
by destabilisation of the proliferative feedback control function h(M). If we compare 
the stem cell model (j3.12|) (with Vq = 1) 

S = b [(1 + X )h (S 1 )S 1 - h (S)S] - e S (4.7) 

with the blast cell model (along the characteristics) 



M 



X)h(M 1 )M 1;1 -h(M)M , (4. 



it is not difficult to sense that modification of the parameters b or A may cause the 
blast maturation to proceed unstably. 
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Figure 6: Proliferatively controlled oscillations due to increased proliferation. Default 
parameter values are used, except that A = 0.6. 



This is what we find if A is increased to 0.6, and the consequent oscillations are 
shown in figure |3 The steady exponential proliferation of blast cells is unstable, 
and this causes oscillations to occur in the maturation profile, and these oscillations 
propagate along the characteristics, as shown in figure 

The oscillations have period equal to the cell cycling time, equal to one in our 
scaled model. A partial understanding of these oscillations is afforded by the obser- 
vation that if h is constant and M is periodic with period 2tt/u, then ()4.8|) admits a 
solution 

M = J2c pq e a ^ +ipuj{t -^\ (4.9) 
p,q 

provided a q satisfies 

a = -A - Ge- a , (4.10) 

where 

A = bh, G = -bh(l + X). (4.11) 

Since A > 0, we have G + A < 0, and it is straightforward to show that there is 
always a single positive root, which can be labelled with q = 0. The others are 
complex (conjugates), and are labelled with increasing frequency as q — ±1, ±2, etc. 
Consideration of these complex roots then shows that for small \G\, Re o q < 0, so that 
the effect of the oscillations dies away as the cells mature; this is what happens in 
figure El However, for larger |G|, Rea q > 0, and the oscillations grow in amplitude as 
the cells mature. This causes M to fluctuate, and thus also h, presumably entraining 
the period of the oscillations to that of the delay. This description is consistent with 
what is seen in figure U\ (see also figure EJ. An approximate criterion for growth of 



B 

M f 
S 
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Figure 7: Two snapshots of the maturation profile for the numerical solution in figure 
El An exponentially growing travelling wave propagates down the maturation axis. 



periodic perturbations with £ is when G < —[A 2 + 



7T 



211/2 



i.e., 



A > 



bh ) 



1/2 



- 1. 



(4.12) 



Differentiation-controlled oscillations 

The final kind of oscillation that we see is induced by the peripheral control of stem 
cell committal through the function Vq(B). These are essentially delay induced oscil- 
lations, where now the delay involved is the maturation time. Because we suppose 
maturation time is large, these are long period oscillations. They can be caused by 
increasing the sensitivity of the peripheral controller, as shown in figure |H1 

To understand the origin of these oscillations, let us suppose that £/ ^> 1, or 
£f ^> r, meaning that the maturation time is significantly longer than the cell cycle 
time, or equivalently that there are a large number of generations in the cell lineage. 
Let us define 

(4.13) 



and the slow time and maturation scales 



We also define /i via 



T = et, X = . 



A = efi, 



(4.14) 



(4.15) 



and suppose that \i = 0(1). Essentially we are revisiting the relaxation oscillation 
analysis of Fowler and Mackey (2002). The partial differential equation for M takes 
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Figure 8: The same graph as in figure [7[ except that a logarithmic scale is used. The 
exponential increase with £ is clearly seen. 




Figure 9: Differentiation-controlled oscillations due to enhancement of the peripheral 
controller function. Default parameters are used, except that v* = v' = 2. 



19 



the form 

dM dM b[h £ , £ Me,e-hM] 

W + dX= e + ^ M *>* 

and expanding in a Taylor series, we have 

d[(l + bh)M] d[(l + bh)M] „ w 

with the boundary condition (taking S = 1) 

M = i/u (S) at X = 0. 



(4.16) 



(4.17) 



(4.18) 



If we suppose h is constant (it is not, but it is not the dependence of h on M 
which causes the oscillations), then the solution of this is 



M = vv [B(T - X)]exp 



jibhX 
1 + bh 



(4.19) 



and the cell efflux at X = 1 (£ = £f) is 

M(l) = i/av [.B(T - 1)], (4.20) 
where the amplification factor a is 

-exp[^y. (4.21) 

Therefore the blood cell conservation law (J3.19)) becomes the delay recruitment model 

JB 



e8 



dT 



uav (Bi 



B. 



(4.22) 



This is a standard delay recruitment equation with a unique steady state. Oscillations 
will occur as a consequence of instability if there are solutions a of (|4.10p . i.e., 



a = -A - Ge~ a , 

with positive real part. The values of A and G are 

A=- G 

so 



ua\v. 



01 



(4.23) 



(4.24) 



The equation ()4.23|) is well understood, see for example Mackey (1978) or Murray 
(2002, pp. 23-26). It is a transcendental equation with an infinite number of complex 
roots which accumulate at the essential singularity at a = — oo. It follows that the 
set of Re a is bounded above. Consequently, there is an instability criterion which 
determines when all the roots a have negative real part, and this is indicated in figure 
HTJl The three curves in the figure are given by G = —A, G = exp[— (1 + A)], and the 
Hopf bifurcation curve G = Gq(A), which is given parametrically by 



.4 



n 

tanfT 



Gq 



sinf2 



(4.25) 
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G 



A 

G = -A 

-2 -1.5 -1 -0.5; 0.5 1 1.5 2 



Figure 10: Stability map for (|4.23jl . The plus and minus signs indicate the sign of 
real values of the growth rate a, when these exist. A Hopf bifurcation occurs as G 
increases through Gq(A), and Go ~ A for large A. 



where Q G [0, 7r]. Since G and A are positive, oscillatory instability occurs precisely 
if G > G (A). Since G ~ A as A — > oo, the instability criterion for large A is simply 
G > A, i.e., 

> 1. (4.26) 
Instability is promoted by increasing \v' \, for example, as indicated in figure El 

5 Conclusions 

In this paper we have studied the onset of oscillations in a model of blood cell pro- 
duction which includes a description of cell cycling and proliferation, and also of 
differentiation and maturation. The model formulation extends the work of previous 
authors by correcting an apparent inconsistency in the description of the primitive 
stem cell population, and also by including the simultaneous control of stem cell pro- 
liferation, stem cell committal, and blast cell proliferation. All three controls can 
cause oscillations for appropriate values of control parameters. 

Previous results concerning stem cell oscillations are reproduced (see figure EJ) but 
these oscillations are harder to obtain when the parameter v is small, and in addition 
they hardly affect the mature blood cell population, without additional destabilisation 
of the blast cell proliferation. The reason for this is that an 0(1) oscillation in the stem 
cell population only causes an 0{y) oscillation in the blast cell committal rate, and 
this amplitude propagates through the differentiating cells. Thus one consequence of 




v\Vq \ exp 



fibh 
1 + bh 
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Figure 11: The effect of switching on all three instability mechanisms. Default pa- 
rameters are used, except that v* = v' = 2, Aq = 0.05, and A = 0.6. 



stem cell paucity is that any instability in the stem cell population itself is hardly 
manifested in the blood cell production. From the point of view of survival and 
control, this is of course a positive result. 

Instability in the proliferation of blast cells due to enhancement of the proliferative 
controller h(M) causes oscillations which propagate down the maturation axis, and 
are amplified as they progress. The result of this is shown in figures El El and |H1 The 
oscillations have a period equal to the cell cycling time. The mechanism for these 
oscillations appears to be a destabilisation of the maturing cell amplification, together 
with a type of resonance which ties the period to the delay. 

The final kind of oscillation is induced by enhanced peripheral control, as seen in 
figure El Stem cell paucity implies that v -C 1, and consequently that > 1, and 
thus that the oscillation period (controlled by the delay is long. This allows an 
approximate reduction of the partial differential delay equation to a simple first order 
differential delay equation, which is simply analysed. In particular, if a threshold 
form of peripheral controller is used, blood cell counts can decrease to zero, as can 
be the case in cyclical neutropenia. 

Finally, and as shown in figure ^2 a combination of all three destabilising mecha- 
nisms can lead to oscillations which operate on both the slow, peripherally controlled 
time scale and the fast, proliferatively controlled one. We consider this observation 
to be a possible explanation of the apparent fact in figure El that both reticulocytes 
and platelets appear to oscillate on a fast as well as a slow time scale. Further study 
of this behaviour requires the extension of this model to accommodate multiple cell 
lineages. 
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